home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
newmat03.lha
/
newmat03
/
newmat3.cxx
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-08
|
11KB
|
381 lines
//$$ newmat3.cxx Matrix get and restore rows and columns
// Copyright (C) 1991: R B Davies and DSIR
#include "include.hxx"
#include "newmat.hxx"
#include "newmatrc.hxx"
//#define REPORT { static ExeCounter ExeCount(__LINE__,3); ExeCount++; }
#define REPORT {}
//#define MONITOR(what,storage,store) \
// { cout << what << " " << storage << " at " << (long)store << "\n"; }
#define MONITOR(what,store,storage) {}
// Control bits codes for GetRow, GetCol, RestoreRow, RestoreCol
//
// LoadOnEntry:
// Load data into MatrixRow or Col dummy array under GetRow or GetCol
// StoreOnExit:
// Restore data to original matrix under RestoreRow or RestoreCol
// IsACopy:
// Set by GetRow/Col: MatrixRow or Col array is a copy
// DirectPart:
// Load or restore only part directly stored; must be set with StoreOnExit
// Still have decide how to handle this with symmetric
// StoreHere:
// used in columns only - store data at supplied storage address, adjusted
// for skip; used for GetCol, NextCol & RestoreCol. No need to fill out
// zeros.
// These will work as a default
// but need to override NextRow for efficiency
void GeneralMatrix::NextRow(MatrixRowCol& mrc)
{
REPORT
if (mrc.cw*StoreOnExit) { REPORT this->RestoreRow(mrc); }
if (mrc.cw*IsACopy)
{
REPORT
real* s = mrc.store + mrc.skip;
MONITOR("Free (NextRow)",mrc.storage,s)
delete [mrc.storage] s;
}
mrc.rowcol++;
if (mrc.rowcol<nrows) { REPORT this->GetRow(mrc); }
else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
}
void GeneralMatrix::NextCol(MatrixRowCol& mrc)
{
REPORT // 423
if (mrc.cw*StoreOnExit) { REPORT this->RestoreCol(mrc); }
if (mrc.cw*IsACopy && (mrc.cw*StoreHere==0))
{
REPORT // not accessed
real* s = mrc.store + mrc.skip;
MONITOR("Free (NextCol)",mrc.storage,s)
delete [mrc.storage] s;
}
mrc.rowcol++;
if (mrc.rowcol<ncols) { REPORT this->GetCol(mrc); }
else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
}
// routines for matrix
void Matrix::GetRow(MatrixRowCol& mrc)
{
REPORT
mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=ncols;
mrc.store=store+mrc.rowcol*ncols;
}
void Matrix::GetCol(MatrixRowCol& mrc)
{
REPORT
mrc.skip=0; mrc.storage=nrows;
if (ncols==1 && (mrc.cw*StoreHere)==0)
{ REPORT mrc.cw-=IsACopy; mrc.store=store; } // not accessed
else
{
mrc.cw+=IsACopy; real* ColCopy;
if ((mrc.cw*StoreHere)==0)
{
REPORT
ColCopy = new real [nrows]; MatrixErrorNoSpace(ColCopy);
MONITOR("Make (MatGetCol)",nrows,ColCopy)
mrc.store = ColCopy;
}
else { REPORT ColCopy = mrc.store; }
if (mrc.cw*LoadOnEntry)
{
REPORT
real* Mstore = store+mrc.rowcol; int i=nrows;
while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
}
}
}
void Matrix::RestoreCol(MatrixRowCol& mrc)
{
// if (mrc.cw*StoreOnExit)
REPORT // 429
if (mrc.cw*IsACopy)
{
REPORT // 426
real* Mstore = store+mrc.rowcol; int i=nrows; real* Cstore = mrc.store;
while (i--) { *Mstore = *Cstore++; Mstore+=ncols; }
}
}
void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); } // 1808
void Matrix::NextCol(MatrixRowCol& mrc)
{
REPORT // 632
if (mrc.cw*StoreOnExit) { REPORT RestoreCol(mrc); }
mrc.rowcol++;
if (mrc.rowcol<ncols)
{
if (mrc.cw*LoadOnEntry)
{
REPORT
real* ColCopy = mrc.store;
real* Mstore = store+mrc.rowcol; int i=nrows;
while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
}
}
else { REPORT mrc.cw -= StoreOnExit; }
}
// routines for diagonal matrix
void DiagonalMatrix::GetRow(MatrixRowCol& mrc)
{
REPORT
mrc.skip=mrc.rowcol; mrc.cw-=IsACopy; mrc.storage=1; mrc.store=store;
}
void DiagonalMatrix::GetCol(MatrixRowCol& mrc)
{
REPORT
mrc.skip=mrc.rowcol; mrc.storage=1;
if (mrc.cw*StoreHere)
{ REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); mrc.cw+=IsACopy; }
else { REPORT mrc.store = store; mrc.cw-=IsACopy; } // not accessed
}
void DiagonalMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
// 800
void DiagonalMatrix::NextCol(MatrixRowCol& mrc)
{
REPORT
if (mrc.cw*StoreHere)
{
if (mrc.cw*StoreOnExit)
{ REPORT *(store+mrc.rowcol)=*(mrc.store+mrc.rowcol); }
mrc.IncrDiag();
if (mrc.cw*LoadOnEntry && mrc.rowcol < ncols)
{ REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); }
}
else { REPORT mrc.IncrDiag(); } // not accessed
}
// routines for upper triangular matrix
void UpperTriangularMatrix::GetRow(MatrixRowCol& mrc)
{
REPORT
int row = mrc.rowcol; mrc.skip=row; mrc.cw-=IsACopy;
mrc.storage=ncols-row; mrc.store=store+(row*(2*ncols-row-1))/2;
}
void UpperTriangularMatrix::GetCol(MatrixRowCol& mrc)
{
REPORT
mrc.skip=0; mrc.cw+=IsACopy; int i=mrc.rowcol+1; mrc.storage=i;
real* ColCopy;
if ((mrc.cw*StoreHere)==0)
{
REPORT // not accessed
ColCopy = new real [i]; MatrixErrorNoSpace(ColCopy);
MONITOR("Make (UT GetCol)",i,ColCopy)
mrc.store = ColCopy;
}
else { REPORT ColCopy = mrc.store; }
if (mrc.cw*LoadOnEntry)
{
REPORT
real* Mstore = store+mrc.rowcol; int j = ncols;
while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
}
}
void UpperTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
{
// if (mrc.cw*StoreOnExit)
{
REPORT
real* Mstore = store+mrc.rowcol; int i=mrc.rowcol+1; int j = ncols;
real* Cstore = mrc.store;
while (i--) { *Mstore = *Cstore++; Mstore += --j; }
}
}
void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
// 722
// routines for lower triangular matrix
void LowerTriangularMatrix::GetRow(MatrixRowCol& mrc)
{
REPORT
int row=mrc.rowcol; mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=row+1;
mrc.store=store+(row*(row+1))/2;
}
void LowerTriangularMatrix::GetCol(MatrixRowCol& mrc)
{
REPORT
int col=mrc.rowcol; mrc.skip=col; mrc.cw+=IsACopy;
int i=nrows-col; mrc.storage=i; real* ColCopy;
if ((mrc.cw*StoreHere)==0)
{
REPORT // not accessed
ColCopy = new real [i]; MatrixErrorNoSpace(ColCopy);
MONITOR("Make (LT GetCol)",i,ColCopy)
mrc.store = ColCopy-col;
}
else { REPORT ColCopy = mrc.store+col; }
if (mrc.cw*LoadOnEntry)
{
REPORT
real* Mstore = store+(col*(col+3))/2;
while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
}
}
void LowerTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
{
// if (mrc.cw*StoreOnExit)
{
REPORT
int col=mrc.rowcol; real* Cstore = mrc.store+col;
real* Mstore = store+(col*(col+3))/2; int i=nrows-col;
while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
}
}
void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
//712
// routines for symmetric matrix
void SymmetricMatrix::GetRow(MatrixRowCol& mrc)
{
REPORT //571
mrc.skip=0; int row=mrc.rowcol;
if (mrc.cw*DirectPart)
{
REPORT
mrc.cw-=IsACopy; mrc.storage=row+1; mrc.store=store+(row*(row+1))/2;
}
else
{
mrc.cw+=IsACopy; mrc.storage=ncols;
real* RowCopy = new real [ncols]; MatrixErrorNoSpace(RowCopy);
MONITOR("Make (SymGetRow)",ncols,RowCopy)
mrc.store = RowCopy;
if (mrc.cw*LoadOnEntry)
{
REPORT // 544
real* Mstore = store+(row*(row+1))/2; int i = row;
while (i--) *RowCopy++ = *Mstore++;
i = ncols-row;
while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
}
}
}
// need to check this out under StoreHere
void SymmetricMatrix::GetCol(MatrixRowCol& mrc)
{
REPORT
mrc.skip=0; int col=mrc.rowcol;
if (mrc.cw*DirectPart)
{
REPORT // not accessed
mrc.cw-=IsACopy; mrc.storage=col+1; mrc.store=store+(col*(col+1))/2;
}
else
{
mrc.cw+=IsACopy; mrc.storage=ncols; real* ColCopy;
if ((mrc.cw*StoreHere)==0)
{
REPORT // not accessed
ColCopy = new real [ncols]; MatrixErrorNoSpace(ColCopy);
MONITOR("Make (SymGetCol)",ncols,ColCopy)
mrc.store = ColCopy;
}
else { REPORT ColCopy = mrc.store; }
if (mrc.cw*LoadOnEntry)
{
REPORT
real* Mstore = store+(col*(col+1))/2; int i = col;
while (i--) *ColCopy++ = *Mstore++;
i = ncols-col;
while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
}
}
}
//void SymmetricMatrix::RestoreRow(int row, real* Rstore)
//{
//// if (cw*IsACopy && cw*StoreOnExit)
// {
// real* Mstore = store+(row*(row+1))/2; int i = row+1;
// while (i--) *Mstore++ = *Rstore++;
// }
//}
//void SymmetricMatrix::RestoreCol(int col, real* Cstore)
//{
//// if (cw*IsACopy && cw*StoreOnExit)
// {
// real* Mstore = store+(col*(col+3))/2;
// int i = nrows-col; int j = col;
// while (i--) { *Mstore = *Cstore++; Mstore+= ++j; }
// }
//}
// routines for row vector
void RowVector::GetCol(MatrixRowCol& mrc)
{
REPORT
mrc.skip=0; mrc.storage=1;
if (mrc.cw >= StoreHere)
{
if (mrc.cw >= LoadOnEntry) { REPORT *(mrc.store) = *(store+mrc.rowcol); }
mrc.cw+=IsACopy;
}
else { REPORT mrc.store = store+mrc.rowcol; mrc.cw-=IsACopy; }
// not accessed
}
void RowVector::NextCol(MatrixRowCol& mrc)
{
REPORT
if (mrc.cw*StoreHere)
{
if (mrc.cw*StoreOnExit) { REPORT *(store+mrc.rowcol)=*(mrc.store); }
// not accessed
mrc.rowcol++;
if (mrc.rowcol < ncols)
{
if (mrc.cw*LoadOnEntry) { REPORT *(mrc.store)=*(store+mrc.rowcol); }
}
else { REPORT mrc.cw -= StoreOnExit; }
}
else { REPORT mrc.rowcol++; mrc.store++; } // not accessed
}
void RowVector::RestoreCol(MatrixRowCol& mrc)
{
REPORT // not accessed
if (mrc.cw>=IsACopy) { REPORT *(store+mrc.rowcol)=*(mrc.store); }
}